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ABSTRACT 

co ' 

We investigate the properties of hybrid gravitational/hydrodynamical 
lO '■ simulations, examining both the numerics and the general physical properties 



CO 



of gravitationally driven, hierarchical collapse in a mixed baryonic/dark matter 
fluid. We demonstrate that, under certain restrictions, such simulations 
converge with increasing resolution to a consistent solution. The dark matter 
achieves convergence provided that the relevant scales dominating nonlinear 
collapse are resolved. If the gas has a minimum temperature (as expected, 
for example, when intergalactic gas is heated by photoionization due to the 
\ ultraviolet background) and the corresponding Jeans mass is resolved, then 

the baryons also converge. However, if there is no minimum baryonic collapse 
mass or if this scale is not resolved, then the baryon results err in a systematic 
fashion. In such resolution is increased the baryon distribution tends 

toward a higher density, more tightly bound state. We attribute this to the 
fact that under hierarchical structure formation on all scales there is always an 
earlier generation of smaller scale collapses, causing shocks which irreversibly 
alter the state of the baryon gas. In a simulation with finite resolution we 
therefore always miss such earlier generation collapses, unless a physical scale is 
introduced below which such structure formation is suppressed in the baryons. 
We also find that the baryon/dark matter ratio follows a characteristic pattern, 
such that collapsed structures possess a baryon enriched core (enriched by 
factors ~ 2 or more over the universal average) which is embedded within a dark 
matter halo, even without accounting for radiative cooling of the gas. The dark 
matter is unaffected by changing the baryon distribution (at least in the dark 
matter dominated case investigated here), allowing hydrodynamics to alter the 
distribution of visible material in the universe from that of the underlying mass. 
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Subject headings: Cosmology: theory — Hydrodynamics — Large scale structure 
of the universe — Methods: numerical 



1. Introduction 



Hydrodynamics is thought to play a key role in the formation of the visible structures 
in the universe, such as bright galaxies and hot intracluster gas. For this reason there is a 
great deal of interest in incorporating hydrodynamical effects into cosmological structure 
formation simulations in order to make direct, quantitative comparisons of such simulations 
to observed data. In addition to gravitation, a cosmological hydrodynamical simulation 
must minimally account for pressure support, shock physics, and radiative cooling, as these 
are the fundamental physical processes thought to play a dominant role in the formation 
of large, bright galaxies (White & Rees 1978). There is already a bewildering array of 
such studies published, including Cen & Ostriker (1992a,b), Katz, Hernquist, & Weinberg 
(1992), Evrard, Summers, & Davis (1994), Navarro & White (1994), and Steinmetz & 
Muller (1994), to name merely a few. In order to appreciate the implications of such 
ambitious studies, it is important that we fully understand both the physical effects of 
hydrodynamics under a cosmological framework and the numerical aspects of the tools 
used for such investigations. Basic questions such as how the baryon to dark matter ratio 
varies in differing structures (galaxies, clusters, and filaments) and exactly how this is 
affected by physical processes such as shock heating, pressure support, or radiative cooling 
remain unclear. It is also difficult to separate real physical effects from numerical artifacts, 
particularly given the current limitations on the resolution which can be achieved. For 
example, in a recent study of X-ray clusters Anninos & Norman (1996) find the observable 
characteristics of a simulated cluster to be quite resolution dependent, with the integrated 
X-ray luminosity varying as L x oc (Ax) _L17 , core radius r c oc (Ax) os , and emission weighted 
temperature Tx oc (Aa;) ' 35 (where Ax is the gridcell size of the simulation). In a study of 
the effects of photoionization on galaxy formation, Weinberg, Hernquist, & Katz (1996) find 
that the complex interaction of numerical effects (such as resolution) with microphysical 
effects (such as radiative cooling and photoionization heating) strongly influences their 
resulting model galaxy population. In this paper we focus on separating physical from 
numerical effects in a series of idealized cosmological hydrodynamical simulations. This 
study is intended to be an exploratory survey of hydrodynamical cosmology, similar in spirit 
to the purely gravitational studies of Melott & Shandarin (1990), Beacom et al. (1991), and 
Little, Weinberg, k Park (1991). 
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We will examine the effects of pressure support and shock heating in a mixed 
baryonic/dark matter fluid undergoing gravitationally driven hierarchical collapse. This 
problem is approached with two broad questions in mind: how stable and reliable is the 
numerical representation of the system, and what can we learn about the physics of such 
collapses? These questions have been investigated for purely gravitational systems in 
studies such as those mentioned above. In those studies numerically it is found that the 
distribution of collisionless matter converges to consistent states so long as the nonlinear 
collapse scale is resolved. Such convergence has not been demonstrated for collisional 
systems, however. It is not clear that hydrodynamical simulations will demonstrate such 
convergence in general, nor if they do that the nonlinear scale is the crucial scale which 
must be resolved. Hydrodynamical processes are dominated by localized interactions on 
small scales, allowing the smallest scales to substantially affect the state of the baryonic 
gas. As an example, consider a collisional fluid undergoing collapse. Presumably such a 
system will undergo shocking near the point of maximal collapse, allowing a large fraction 
of the kinetic energy of the gas to be converted to thermal energy. In a simple case such as 
a single plane- wave perturbation (the Zel'dovich pancake collapse), the obvious scale which 
must be resolved is the scale of the shock which forms around the caustic. However, in a 
hierarchical structure formation scenario there is a hierarchy of collapse scales, and for any 
given resolution limit there is always a smaller scale which will undergo nonlinear collapse. 
The subsequent evolution of the gas could well depend upon how well such small scale 
interactions are resolved, and changes in the density and temperature of gas on small scales 
could in turn influence how it behaves on larger scales (especially if cooling is important). 

In this paper we examine a series of idealized experiments, evolving a mixed fluid 
of baryons and collisionless dark matter (dark matter dominated by mass), coupled 
gravitationally in a flat, Einstein-de Sitter cosmology. The mass is seeded with Gaussian 
distributed initial density perturbations with a power-law initial power spectrum. We 
perform a number of simulations, varying the resolution, the initial cutoff in the density 
perturbation spectrum, and the minimum allowed temperature for the baryons. Enforcing 
a minimum temperature for the baryons implies there will be a minimal level of pressure 
support, and therefore a minimum collapse scale (the Jeans mass), below which the baryons 
are pressure supported against collapse. From the numerical point of view, performing 
a number of simulations with identical initial physical conditions but varying resolution 
allows us to unambiguously identify resolution effects. By enforcing a Jeans mass for the 
baryons we introduce an intrinsic mass scale to the problem, which may or may not be 
resolved in any individual experiment. The hope is that even if the gas dynamical results do 
not converge with increasing resolution in the most general case, the system will converge if 
the fundamental Jeans mass is resolved. 
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The effects of the presence (or absence) of a baryonic Jeans mass also raises interesting 
physical questions. Though we simply impose arbitrary minima for the baryon temperatures 
here, processes such as photoionization enforce minimum temperatures in the real universe 
by injecting thermal energy into intergalactic gas. The Gunn-Peterson test indicates that 
the intergalactic medium is highly ionized (and therefore at temperatures T ^ 10 4 K) out 
to at least z <^ 5. Shapiro, Giroux, & Babul (1994) discuss these issues for the intergalactic 
medium. The dark matter, however, is not directly influenced by this minimal pressure 
support in the baryons, and therefore is capable of collapsing on arbitrarily small scales. 
Pressure support provides a mechanism to separate the two species, and since the dark 
matter dominates the mass density it can create substantial gravitational perturbations on 
scales below the Jeans mass. While there are many studies of specific cosmological models 
with detailed microphysical assumptions, the general problem of the evolution of pressure 
supported baryons in the presence of nonlinear dark matter starting from Gaussian initial 
conditions has not been investigated in a systematic fashion. 

This paper is organized as follows. In §[| we discuss the particulars of how the 
simulations are constructed and performed. In ^ we characterize the numerical effects 
we find in these simulations, and in §|4] we discuss our findings about the physics of this 
problem. Finally, §|5| summarizes the major results of this investigation. 

2. The Simulations 

A survey such as this optimally requires a variety of simulations in order to adequately 
explore the range of possible resolutions and input physics. Unfortunately, hydrodynamical 
cosmological simulations are generally quite computationally expensive, and therefore in 
order to run a sufficiently broad number of experiments we restrict this study to 2-D 
simulations. There are two primary advantages to working in 2-D rather than 3-D. First, 
parameter space can be more thoroughly explored, since the computational cost per 
simulation is greatly reduced and a larger number of simulations can be performed. Second, 
working in 2-D enables us to perform much higher resolution simulations than are feasible 
in 3-D. While the real universe is 3-D and we must therefore be cautious about making 
specific quantitative predictions based on this work, 2-D experiments can be used to yield 
valuable qualitative insights into the behaviour of these systems. For similar reasons Melott 
& Shandarin (1990) and Beacom et al. (1991) also utilize 2-D simulations in their studies 
of purely gravitational dynamics. 

The 2-D simulations presented here can be interpreted as a slice through an infinite 
3-D simulation (periodic in (x,y) and infinite in z). The particles interact as parallel rods 
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of infinite length, obeying a gravitational force law of the form F grav oc 1/r. The numerical 
technique used for all simulations is SPH (Smoothed Particle Hydrodynamics) for the 
hydrodynamics and PM (Particle-Mesh) for the gravitation. The code and technique are 
described and tested in Owen et al. (1996), so we will not go into much detail here. We 
do note, however, that while our code implements ASPH (Adaptive Smoothed Particle 
Hydrodynamics) as described in our initial methods paper, we are not using the tensor 
smoothing kernel of ASPH for this investigation, but rather simple SPH. The results should 
be insensitive to such subtle technique choices since the goal is to compare simulation to 
simulation, so we employ simple SPH in order to separate our findings from questions of 
technique. 

All simulations are performed under a flat, Einstein-de Sitter cosmology, with 10% 
baryons by mass (f2 ba ry = 0.1,fidm = 0.9, A = 0). Thus the mass density is dominated by 
the collisionless dark matter, which is linked gravitationally to the collisional baryons. The 
baryon and dark matter particles are initialized on the same perturbed grid, with equal 
numbers of both species. Therefore, initially all baryons exactly overlie the dark matter 
particles, and only hydrodynamical effects can separate the two species. The baryon/dark 
matter mass ratio is set by varying the particle mass associated with each species. The 
initial density perturbation spectrum is taken to be a power-law P{k) = (\5p(k)/p\ 2 ) oc k n 
up to a cutoff frequency k c . Note that since these are 2-D simulations, for integrals over the 
power spectrum this is equivalent in the 3-D to a power spectrum of index n — 1. In this 
paper we adopt a "flat" (n = 0) 2-D spectrum 



where A norm normalizes the power-spectrum. Note that using a flat cosmology and 
power-law initial conditions implies these simulations are scale-free, and should evolve 
self-similarly in time. We can choose to assign specific scales to the simulations in order 
to convert the scale-free quantities to physical units. All simulations are halted after 60 
expansion factors, at which point the nonlinear scale (the scale on which Sp/p ~ 1) is 
roughly 1/8 of the box size. 

The Jeans length is the scale at which pressure support makes the gas stable against 
the growth of linear fluctuations due to self-gravitation - the Jeans mass is the amount of 
mass contained within a sphere of diameter the Jeans length. The Jeans length Aj and 
mass Mj are defined by the well known formula (Binney & Tremaine 1987) 




(1) 
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where p is the mass density and c s the sound speed. The baryons are treated as an ideal 
gas obeying an equation of state of the form P = (7 — l)up, where P is the pressure and u 
is the specific thermal energy. Enforcing a minimum specific thermal energy (and therefore 
temperature) in the gas forces a minimum in the sound speed c 2 = jP/ p = 7(7 — l)u, which 
therefore implies we have a minimum Jeans mass through equation (^|). Note that p is the 
total mass density (baryons and dark matter), since it is the total gravitating mass which 
counts, and therefore Mj as expressed in equation (^) represents the total mass contained 
within r < Xj/2. If we want the total baryon mass contained within this radius, we must 
multiply Mj by fi b ar y /^- 

It is also important to understand how the mass resolution is set for the baryons 
by the SPH technique. This is not simply given by the baryon particle mass, since SPH 
interpolation is a smoothing process typically extending over spatial scales of several 
interparticle spacings. In general the mass resolution for the hydrodynamic calculations can 
be estimated as the amount of mass enclosed by a typical SPH interpolation volume. If the 
SPH smoothing scale is given by h and the SPH sampling extends for 77 smoothing scales, 
then the mass resolution Mr is given by 

M R = ^(rjhfp. (4) 

This is probably something of an overestimate, since the weight for each radial shell in this 
interpolation volume (given by the SPH sampling kernel W) falls off smoothly towards 
r = 7]h, but given the other uncertainties in this quantity equation (f|) seems a reasonable 
estimate. Note that the resolution limit for the SPH formalism is best expressed in terms of 
a mass limit, appropriate for SPH's Lagrangian nature. For this reason we choose to express 
the Jeans limit in terms of the Jeans mass (eq. ||) throughout this work, as the Jeans 
limit can be equally expressed in terms of a spatial or a mass scale. In N-body work it is 
common to express the mass resolution of an experiment in units of numbers of particles. In 
our simulations we use a bi-cubic spline kernel which extends to 77 = 2 smoothing lengths, 
and initialize the smoothing scales such that the smoothing scale h extends for two particle 
spacings. We therefore have a mass resolution in 2-D of roughly 50 particles, or equivalently 
in 3-D roughly 260 particles. 

We perform simulations both with and without a minimum temperature (giving Jeans 
masses Mj = 0, Mj > 0), at three different resolutions (N = iV bary = iV dm = 64 2 , 128 2 , 
and 256 2 ), and for three different cutoffs in the initial perturbation spectrum (k c = 32, 64, 
and 128). The initial density perturbations are initialized as Gaussian distributed with 
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random phases and amplitudes, but in such a manner that all simulations have identical 
phases and amplitudes up to the imposed cutoff frequency k c . The cutoff frequencies 
are the subset of k c G (32, 64, 128) up to the Nyquist frequency for each resolution 
fc Nyq = N l ' 2 /2, so for each resolution we have k c (N = 64 2 ) = 32, k c {N = 128 2 ) G (32,64), 
and k c (N = 256 2 ) G (32, 64, 128). For each value of the minimum temperature we therefore 
have a grid of simulations which either have the same input physics at differing resolutions 
(i.e., k c = 32 for N G [64 2 , 128 2 , 256 2 ]), or varying input physics at fixed resolution (i.e., 
N = 256 2 for k c G [32,64,128]). This allows us to isolate and study both numerical 
and physical effects during the evolution of these simulations. In total we discuss twelve 
simulations. 

For the simulations with a minimum temperature, there is an ambiguity in assigning 
a global Jeans mass with that temperature. The density in equation (^) is formally the 
local mass density, and therefore the Jeans mass is in fact position dependent through p(r). 
Throughout this work we will refer to the Jeans mass at any given expansion factor as 
the Jeans mass defined using a fixed minimum temperature and the average background 
density, making this mass scale a function of time only. This is equivalent to taking the 
zeroth order estimate of Mj, giving us a well defined characteristic mass scale. In terms of 
this background density, Figure [I] shows the baryon Jeans mass (in units of the resolved 
mass via equation [[|) as a function of expansion. Note that for a given simulation Mr 
remains fixed, and it is the Jeans mass which grows as Mj oc p~ l l 2 oc a 3 / 2 . It is apparent 
that the N = 256 2 simulations resolve the Jeans mass throughout most of the evolution, 
the N = 128 2 simulations resolve Mj by a/at ~ 15, and the N = 64 2 simulation does not 
approach Mj/Mr ~ 1 until the end of our simulations at ajai ~ 60. The specific value 
of T min used in this investigation is chosen to yield this behaviour. We discuss physically 
motivated values for this minimum temperature in §|5|. 



3. Numerical Resolution and the Jeans Mass 
3.1. Dark Matter 

We will begin by examining the dark matter distribution, as this is a problem which 
has been examined previously. Figures ^|a and b show images of the dark matter overdensity 
Pdm/pdm for the Mj = simulations. In order to fairly compare with equivalent images 
of the SPH baryon densities, the dark matter information is generated by assigning a 
pseudo-SPH smoothing scale to each dark matter particle, such that it samples roughly the 
same number of neighboring dark matter particles as the SPH smoothing scale samples 
in the baryons. We then use the normal SPH summation method to assign dark matter 
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densities, which are used to generate these images. The panels in the figure are arranged 
with increasing simulation resolution N along rows, and increasing cutoff frequency k c down 
columns. The diagonal panels represent each resolution initialized at its Nyquist frequency 
for P(k). Note that the physics of the problem is constant along rows, and numerics is 
constant along columns. If resolution were unimportant, the results along rows should be 
identical. Likewise, since the numerics is held constant along columns, only physical effects 
can alter the results in this direction. 

Comparing the dark matter densities along the rows of Figure 0a, it is clear that 
the structure becomes progressively more clearly defined as the resolution increases. This 
is to be expected, since the higher resolution simulations can resolve progressively more 
collapsed/higher density structures. The question is whether or not the underlying particle 
distribution is systematically changing with resolution. In other words, do the simulations 
converge to the same particle distribution on the scales which are resolved? Figure 0b 
shows this same set of dark matter overdensities for the Mj = simulations, only this 
time each simulation is degraded to an equivalent N = 64 2 resolution and resampled. 
This is accomplished by selecting every nth node from the higher resolution simulations, 
throwing away the rest and suitably modifying the masses and smoothing scales of the 
selected particles. Note that now the dark matter distributions look indistinguishable for 
the different resolution experiments, at least qualitatively. This similarity implies that the 
high frequency small scale structure has minimal effect on the larger scales resolved in this 
figure. Looking down the columns of Figure 0a it is clear that increasing k c does in fact 
alter the dark matter particle distribution, such that the large scale, smooth filaments are 
progressively broken up into smaller clumps aligned with the overall filamentary structure. 
These differences are lost in the low-res results of Figure 0b, implying that these subtle 
changes do not significantly affect the large scale distribution of the dark matter. 

In Figures |3|a and b we show the mass distribution functions for the dark matter 
overdensity /(pdm/pdm)- Figure [5]a includes all particles from each simulation (as in Figure 
0a), while Figure ||b is calculated for each simulation degraded to equivalent N = 64 2 
resolutions (comparable to Figure |||b). The panels are arranged as in Figure 0, with Mj = 
and Mj > overplotted as different line types. It is clear that the varying Jeans mass in 
the baryons has negligible effect on the dark matter, a point we will return to in §[|. The 
full resolution results of Figure 0a show a clear trend for a larger fraction of the mass to 
lie at higher densities with increasing resolution. There is also a similar though weaker 
trend with increasing k c . However, examining the resampled results of Figure ^]b it appears 
that the results of all simulations converge, bearing out the visual impressions of Figures 
^|a and b. For the dark matter, with increasing resolution more information is gained 
about the highest density/most collapsed fraction of the mass, but so long as the pertinent 
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nonlinear scales are resolved the results converge. The underlying particle distribution does 
not depend upon the numerical resolution, similarly to the results discussed in Little et al. 
(1991). 

3.2. Baryons 

We now turn our attention to the baryon distribution. Figures f|a, b, c, and d show 
images of the baryon overdensity for Mj = and Mj > at expansions ajai = 30 and 
a/cii = 60. There is a pronounced trend for the collapsed filaments and clumps to become 
progressively more strongly defined as the simulation resolution improves - even more so 
than we see in the dark matter. The tendency to break up filaments into small scale clumps 
with increasing k c is also clearly evident for the Mj = case. Additionally, the presence of 
a nonzero Jeans mass visibly influences the baryon density distribution in Figures [|c and 
d. This is particularly evident in the high resolution iV = 256 2 column, where the increased 
pressure support creates a "puffier" distribution, wiping out the smallest scale structures 
in the baryons. Recall from Figure |1] that we naively expect the presence of the pressure 
support for Mj > to affect both N = 128 2 and N = 256 2 at a/a, t = 30, but not N = 64 2 . 
Comparing the results of Figures [|a and c, we indeed see this trend. By a/a* = 60, the 
effects of the Jeans mass are clearly evident (comparing Figures [§b and d) for N = 128 2 
and N = 256 2 , though N = 64 2 still appears relatively unaffected. 

Figures |5|a and b show images of the baryon densities for the k c = 32 simulations, but in 
this case resampled to iV = 64 2 resolutions analogous to Figure |2|b. At expansion a/a,i = 30 
(Figure ^]a), we see that for Mj = the baryons appear to be systematically more tightly 
collapsed with increasing simulation resolution, even though they have all been resampled 
to the same sampling resolution to produce this image. This supports the view that the 
baryon distribution is fundamentally changing with increasing simulation resolution, in 
contrast with the dark matter. The N = 128 2 and N = 256 2 Mj > simulations, however, 
demonstrate very similar baryon density images, though iV = 64 2 still appears different at 
a/cii = 30. At a/cii = 60 (Figure |5|b) we again see for Mj = a clear trend with simulation 
resolution, while the Mj > runs look remarkably similar to one another. 

Figures |5|a and b show the full resolution mass distribution functions of the baryon 
overdensities /(pbary/pbary) for all simulations at a/a« = 30 and a/ai = 60, respectively. 
The Mj = functions show a strong trend to transfer mass from low to high densities 
with increasing resolution, and a similar though weaker trend with k c . However, even at 
full resolution the Mj > simulations show very similar density distributions once Mj is 
resolved. The Mj > simulations also appear to be relatively insensitive to k c , suggesting 
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that the increased small scale power is being wiped out by the pressure support. Figures 
^|c and d show these same baryon density distribution functions, only for all simulations 
degraded to iV = 64 2 resolutions. These bear out our previous observations. In the case with 
no Jeans mass, there is no sign of convergence in the baryon distribution as the resolution 
is increased. However, when a Jeans mass is present, then the baryon distributions do 
converge once the Jeans mass is resolved. 

Figure [7| presents a more quantitative way to measure this convergence problem. In 
this figure we calculate the Kolmogorov-Smirnov statistic -D(pbary), comparing the baryon 
density distribution for each simulation to the others at the same expansion and Jeans mass. 
We do not expect these simulations to exactly reproduce one another, and therefore there is 
little point in assigning significance to the exact quantitative value of D. However, the K-S 
statistic does provide objective measures of how similar or dissimilar these distributions 
are, and therefore we might expect to learn something by comparing their relative values. 
Comparing the upper panels of Figure [7] we can see that at a/ai = 30 the N = 128 2 and 
iV = 256 2 simulations are more similar for Mj > than for Mj = 0, while the N = 64 2 
simulation remains relatively distinct in both cases. At a/ai = 60, however, we can see that 
for Mj > all the simulations appear comparable, while for Mj = they remain distinct 
for the different resolutions. 

We therefore have a subtly different picture for the numerical behaviour of the dark 
matter and baryons. The critical resolution scale for the dark matter is the scale of 
nonlinearity. So long as this scale is resolved, the dark matter distribution can be expected 
to converge to a consistent state on resolved scales. Unfortunately, the distribution and 
state of the baryonic particles appears in general to be sensitive to the numerical resolution. 
However, it is possible and physically plausible to define a fundamental collapse scale in 
the form of the Jeans mass for the baryons, below which baryonic structure formation is 
suppressed. This scale can now be treated as the critical baryonic resolution scale, and we 
do find that once this threshold is reached the baryon distribution will reliably converge as 
well. 



4. Hydrodynamics and the Baryon Distribution 

4.1. Shocks and Temperatures 

The results of the previous section indicate that hydro dynamical interactions on small 
scales can significantly alter the the final state of the baryons in ways which propagate 
upward and affect larger scales. The tendency for a simulation with a given finite resolution 
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is to underestimate the "true" fraction of high density, collapsed baryons. A likely cause 
for this trend is the presence of small scale, unresolved shocks in the baryon gas. Because 
shocks provide a mechanism for transferring the gas's kinetic energy to thermal energy, it 
is reasonable to expect that the fashion and degree to which the baryons collapse will be 
dependent upon when and how strongly they undergo shocks. In this section we investigate 
the thermal state of the baryons, with the goal of understanding the pattern and importance 
of shocking in the gas. 

In the top row of Figure [8] we show the 2-D mass distribution function of the baryons in 
terms of their overdensity and temperature /(pbary/Pbary, T) for the Mj = simulations at 
a/cii = 60. The various resolutions share some gross properties in the p — T plane. The low 
density gas tends for the most part to be cool, though there is a tail of low density material 
with temperatures up to T ^ 10 4 K. The high density gas is at relatively high temperatures, 
with most of the material near T ~ 10 6 K. However, there is a notable trend for the highest 
density material to be somewhat cooler with increasing simulation resolution. This effect is 
similar to the behaviour seen in simple 1-D collapse such as the Zel'dovich pancake (Shapiro 
& Struck-Marcell 1985). The highest density gas is the fraction which collapses earliest, 
when the background density is highest. Such gas is placed on a lower adiabat than gas 
which falls in at later times, and thus remains cooler. In our case this means that since 
higher resolution simulations can resolve higher density clumps (which therefore form at 
earlier times), we should tend to see the temperature of the highest density material fall 
with increasing resolution. 

The high temperature gas is heated by shocks as it falls into the dark matter dominated 
potential wells. In order to isolate shock heating from simple adiabatic compression heating, 
we calculate the distribution of the temperature in units of the adiabatic temperature T ad , 
given by 

T ad = T (^Y\ (5) 

T ad represents the temperature the gas would be at if it were only heated through simple 
PdV work. Since the only non-adiabatic process we allow is shock heating, only gas which 
has undergone shocking should be at T/T a d > 1. In the bottom row of Figure |8] we calculate 
the distribution /(pbary/Pbary> T/Tad) f° r the Mj = simulations at a/aj = 60. The high 
density fraction of the gas is clearly strongly shocked in all cases, with T/T a( j ~ 10 7 — 10 9 . 
There is a clear trend for T/T ac j in the high density gas to fall with resolution, indicating 
that the highest density fraction of the gas is less strongly shocked as the resolution 
increases. Though we do not show the results at fixed resolution and increasing k c here, 
there are also subtle trends evident with k c in both the p — T and p — T /T & a planes. 
Generally the temperature/shocking distribution of moderately overdense material grows 
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wider with increasing k c . 

It appears that shocks are indeed the key physical mechanism distinguishing the 
different resolution experiments. We find that in general most of the baryonic material is 
processed through shocks at some point. We note a general pattern in which the highest 
density gas in low resolution experiments is characteristically more strongly shocked than 
the highest density gas in higher resolution experiments. The physical inference of these 
trends is that the larger the region which collapses, the stronger the resulting shock. The 
underlying physical mechanism for this property is easily understood. Since the highest 
density material represents the gas which collapses earliest, this is also the gas which falls 
into the shallowest potential wells. As the structure continues to grow, these potential wells 
deepen. Gas which infalls at later times therefore picks up more kinetic energy, which in 
turn leads to stronger shocking and higher temperatures. Once shocking occurs, the state 
of the baryon gas is discontinuously and irreversibly altered. In order to properly represent 
the physical state of the gas, a simulation must resolve the smallest scales on which shocks 
are occuring. This is why enforcing a Jeans mass allows convergence, since establishing a 
minimum Jeans mass implies there is a minimum scale on which baryonic structures can 
form, forcing a minimum scale for shocking. 



4.2. Comparing the Baryon & Dark Matter Distributions 

One of the most fundamental questions we can address is how the distributions of 
dark matter and baryons compare to one another. Comparing the dark matter and baryon 
density fields for the Mj = case in Figures ^a and |]b, there is a distinct impression that 
the baryons tend to be more tightly clustered than the dark matter on all collapsed scales. 
The situation is a bit more complex for the Mj > case in Figure f|d. Comparing the 
(N = 256 2 , k c = 128) distributions, it is evident that the Mj > baryons show a more 
diffuse structure than that of the Mj = case, to the point that some of the smallest scale 
structures are entirely suppressed. Bear in mind that the dark matter evolves essentially 
independently of the baryons in this dark matter dominated case, so the small scale 
structures still form in the overall mass distribution - the baryons are simply excluded from 
them. The large scale structures such as the filaments and the largest knots are still quite 
prominent in the Mj > baryon distribution, just as for the Mj = case. These patterns 
suggest that the baryons are generically more clustered than the dark matter, down to the 
scale set by the Jeans mass. At this scale and lower, the dark matter continues to form 
collapsed structures, whereas the baryons are held out of these structures by the pressure 
support enforced by the minimum temperature. 
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In Figure |9] we calculate the baryon to dark matter number density ratio as 
a function of baryonic overdensity. The baryon to dark matter ratio is defined as 
rzbary/^dm = ^dmPbary/^baryPdm, so that n haxy / n Ara > 1 corresponds to baryon enrichment, 
while nbary/^dm < 1 implies baryon depletion. There is a clear trend for the highest density 
material to be baryon enriched, implying that the cores of the most collapsed structures are 
relatively enriched in baryons compared with the universal average. This trend persists even 
in the Mj > simulations, though it is not as pronounced as in the Mj = case. There is 
no evidence that a significant fraction of the baryons exist in regions which are dark matter 
enhanced. In all simulations underdense material appears to lie near the universal average 
mixture n^y/rid m ~ 1. We also note a trend with resolution, such that the higher the 
resolution of the simulation, the greater the baryon enrichment found in overdense regions. 

A simple physical picture can account for these trends. So long as the density evolution 
is in the linear regime (Sp/p <C 1), the dark matter and baryons evolve together, remaining 
at the universal mix of nbary/^dm ~ 1- During this linear phase the pressure support 
(barring any imposed minimum pressure) is orders of magnitude less important than the 
gravitational term, so the baryon/dark matter fluid evolves as a pressureless gas. Once 
nonlinear collapse begins {Sp/p ^ 1), the baryons rapidly fall inward with the dark matter 
until they collide near the potential minimum. At this point the baryon gas shocks, 
converting the majority of its kinetic energy into thermal energy, and it stops, forming a 
hot pressure supported gas at the bottom of the potential well. In the case with a minimum 
pressure support, the collapse proceeds until the pressure term (due to the increase in 
density) builds sufficiently to impede the baryons infall, at which point the baryons slow, 
separate from the infall, and shock. In either case the dark matter forms a more diffuse 
structure supported by velocity dispersion. This process leads to the generic patterns noted 
above: on scales below which the collapse has become nonlinear, the baryons tend to be 
characteristically more clustered than the dark matter, at least down to the minimal point 
set by the Jeans scale. In either case the critical factor determining exactly when the 
baryons separate from the general inflow is the point at which shocking sets in. We also 
know from the numerical observations that this process is resolution dependent, and in fact 
the baryon enrichments we see for the Mj = case in Figure ^ must represent lower limits 
to the "true" baryon enrichment. The enrichments noted for the Mj > simulations should 
be reliable, to the extent that the specific minimum temperature chosen is reasonable. 

It is somewhat puzzling to note that our measured positive biasing of the baryons in 
collapsed structures is at odds with previously published results. In a study of the cluster 
formation under the standard Q = 1 Cold Dark Matter (CDM) model, Evrard (1990) finds 
that while outside of the cluster environment the baryons and dark matter simply track the 
universal average mix, the baryon fraction within the cluster is in fact somewhat lowered. 
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Kang et al. 1994 examine a larger volume of an = 1 CDM cosmology, and find that not 
only are the overdense regions baryon depleted, but that their underdense, void structures 
are baryon enriched. There are several possible explanations for this disagreement. One 
possibility is that this represents a geometric effect, in that our experiments are 2-D, 
while these other studies employ fully 3-D simulations. In our simulations, the "filaments" 
actually represent walls, and the most collapsed knots are best interpreted as cross-sections 
through tubular filaments. The processes of collapsing to a plane, a line, or a point are 
certainly different processes, and the isotropy of pressure support makes these structures 
progressively more difficult to form. In a 1-D planar collapse, for instance, it is well known 
that the central collapse plane will be baryon enriched, while the question of whether or 
not a cluster is baryon enriched or depleted is still hotly debated. We see some evidence 
for this effect in Figure |9|. Looking particularly at the upper dashed lines in this figure 
(representing the baryon enrichment at which 90% of the mass at that overdensity lies 
below) we note our most extreme enrichments occur at moderate overdensities, roughly 
in the range Pbary/Pbary ~ 10 1 — 10 2 . This extremely baryon enriched material represents 
the "filaments" in our simulations (walls in 3-D). It is also possible that resolution effects 
play a role here. As pointed out previously, we find a strong resolution dependence, such 
that finite resolution tends to underestimate the fraction of high density, collapsed baryonic 
material. Evrard (1990) uses 16 3 SPH nodes to represent his baryon component, which for 
the scale of his box is equivalent to our lowest resolution simulations. Kang et al. (1994) 
use an entirely different technique to simulate the hydrodynamics, which relies upon a fixed 
grid to represent the baryons. This limits their spatial resolution so that typical clusters are 
only a handful of cells across. It is also important to compare these quantities in the same 
manner. In Figure |9| we calculate the baryon to dark matter mixture in a manner which 
follows the baryon mass, since we sample at the positions of the baryon particles. This 
naturally gives the greatest weight to the most prominent baryonic structures. Kang et al. 
(1994) calculate this distribution in a manner which is volume weighted, which will tend to 
give the greatest weight to underdense, void like regions. Since the baryon fraction appears 
to be a function of environment, these differences can be significant. Without further study, 
it is difficult to know the true reason for this discrepancy, or how the actual baryon/dark 
matter ratio should evolve. 



5. Discussion 

The results of this investigation can be broken into two broad categories: what is 
revealed about the physics of hierarchical collapse in a mixed baryonic/dark matter fluid, 
and what is learned about the numerics of simulations of this process. We find that the dark 
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matter converges to a consistent state on resolved scales, so long as the nonlinear collapse 
scales are well resolved. Increasing the resolution of the experiment does not fundamentally 
alter the dark matter distribution, but simply yields more detailed information about the 
small scale collapsed structures. This is in agreement with previous, purely collisionless 
studies, though we demonstrate this here including a collisional component. 

The numerical story is quite different for the collisional baryonic gas. We find 
that in the case where we do not impose a fundamental physical resolution scale in the 
baryons, the simulation results do not converge with increasing resolution. Rather, as the 
numerical resolution of the experiment is increased, the collapsed fraction of the baryons 
is systematically altered toward a higher density, more tightly bound, and less strongly 
shocked state. The physical reason for this behaviour is the presence of shocks, which 
allows the evolution on small scales to affect the overall state of the baryonic mass. With 
improving resolution the simulation is able to resolve the collapse of smaller structures at 
earlier times. The smaller scale (and therefore earlier) the resolved collapse, the weaker the 
resulting shock is found to be. This effect is most obvious in Figure where there is a 
systematic trend of higher density/more weakly shocked material with increasing resolution. 

The fact that the dark matter converges in general with resolution, while the baryons 
do not, highlights a fundamental difference in the physics of these two species. While 
both dark matter and baryon fluids react to the global and local gravitational potential, 
the baryons are additionally subject to purely local hydrodynamical phenomena - most 
prominently shocking in this case. Once strong shocking sets in these hydrodynamical 
effects can rise to rival the gravitational force on the baryonic fluid, allowing the baryons 
to be strongly influenced by interactions on small scales in ways which the dark matter 
is not. This implies that such small scale interactions can be just as important as the 
large scale forces in determining the final state of the baryons. In other words, for the 
dark matter there is no back reaction from small to large scales, whereas the baryons 
are strongly influenced by interactions on small scales. In the coupling of these physical 
processes, gravitation dominates the large scale structure, but hydrodynamics affects the 
local arrangements and characteristics of the baryonic gas. If we want the quantitative 
results of such studies to be reliable, we must have reason to believe that the localized 
hydrodynamical processes are adequately resolved. 

This gloomy picture is alleviated by an important physical effect: the Jeans mass. 
Introducing a minimum temperature (and therefore pressure support) into the baryons 
creates a fundamental length/mass scale, below which the baryons are supported by pressure 
against any further collapse or structure formation. We find that once we introduce such 
a minimal scale into the baryonic component, the simulation results converge as this scale 
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is resolved. This convergence holds even though the dark matter component continues to 
form structures below the baryon Jeans scale. Although the Jeans scale is dependent upon 
the local density, we find that the global Jeans scale defined using the background density is 
adequate to define the critical resolution necessary for the hydrodynamics to converge. This 
therefore describes an additional resolution scale necessary for hydrodynamical simulations 
to meet, much as the nonlinear mass scale represents the crucial resolution necessary for 
purely gravitational systems. Furthermore, our experiments indicate that equation @ is 
a reasonable estimate of an SPH simulation's true mass resolution, since we find that the 
threshold Mr <; Mj marks the point at which convergence is achieved. 

In these experiments we have tested the effects of the Jeans mass in an idealized 
framework by simply imposing an arbitrary minimum temperature into our system, but 
there is reason to believe that such minimum temperatures should exist in the real universe. 
Based upon observations such as the Gunn- Peterson test (Gunn & Peterson 1965), it is 
known that the IGM is highly ionized out to redshifts z <^ 5, which implies a minimum 
temperature for the IGM of at least T ^ 10 4 . Assuming an Einstein-de Sitter cosmology, a 
minimum temperature of T ~ 10 4 requires a minimum spatial resolution (via eq. @) 

A^ 0.777 (l + ^(A)-" 2 (_Z_) ,/2 ft - lM pc, (6) 
which equates to a baryon mass resolution of (eq. 0) 

/ n \ -3/2 / rp .3/2 

M R & fi bary Mj ~ 6.82 x 10 10 fi bary (1 + z)-'" 2 h- 1 M Q . (7) 

This limit can also be expressed in terms of a minimum circular velocity, which has the 
advantage of being independent of redshift. The minimum circular velocity can found as a 
function of the minimum temperature by relating the kinetic energy necessary for dynamical 
support to the internal energy for equivalent pressure support (Thoul & Weinberg 1996), 
yielding 

(2kT\ 1/2 ( a \-V2 / T \V\ , 

~ 16 - 6 (^) Go*) km/sec - (8) 

In our Mj > simulations if we choose to call the scale at which RMS mass fluctuation is 
AM/M ~ 0.5 to be 8 h^ 1 Mpc at the final expansion, then our box scale is L = 64/i~ 1 Mpc 
and the minimum temperature corresponds to T min ~ 10 6 K. While there are some 
suggestions that the intergalactic medium could be heated to temperatures as hot as 10 6 K 
(through mechanisms such as large scale shocks of the IGM), clearly these simulations 
do not meet our criteria if we wish to consider photoionization as setting the minimum 
temperature. It is also not clear that the current generation of large-scale hydrodynamical 
cosmological simulations meet this criterion, but it should be achievable. 
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It is still unclear whether or not in the case with no minimum temperature imposed the 
baryon distribution will eventually converge. It is well known that in a purely gravitational 
system, as structure builds and smaller dark matter groups merge into larger structures, the 
dark matter "forgets" about the earlier small scale collapses as such small structures are 
incorporated into larger halos and disrupted. This is why the dark matter results converge 
once the nonlinear mass scale is resolved. While it is evident from studies such as this that 
the baryons maintain a longer memory of their previous encounters, it seems likely that 
as the baryon gas is progressively processed through larger scale and stronger shocks, at 
some point the previous evolution should become unimportant. At exactly what level this 
transition is reached remains uncertain, however, as we see no evidence for such convergence 
here. 

Radiative cooling must be accounted for in order to model processes such as galaxy 
formation, and the inclusion of radiative cooling can only exacerbate the non-convergence 
problems we find here. The amount of energy per unit mass dissipated by radiative cooling is 
proportional to the density, and we have already noted that the trend with finite resolution 
is to underestimate the local gas density and overestimate the temperature. Given these 
tendencies, it is not difficult to envision problems for finite resolution simulations which 
will tend to underestimate the effectiveness of radiative cooling in lowering the temperature 
(and therefore pressure support) of the shocked gas. This could lead to perhaps drastic 
underestimates of the fraction of cold, collapsed baryons for a given system, and therefore 
strongly influence the inferred galaxy formation. Evrard et al. (1994) note this effect when 
comparing their high and low resolution 3-D simulations. They find that altering their 
linear resolution by a factor of two (and therefore the mass resolution by a factor of eight) 
changes the measured total amount of cold collapsed baryons by a factor of ~ 3. They 
attribute this change to just the sort of problems we discuss here. Weinberg, Hernquist, & 
Katz (1996) report similar findings and interpretation for simulations with a photoionizing 
background. It therefore seems likely it is all the more important to resolve the minimum 
mass scale set by the minimum temperature in systems with radiative cooling. 

We find that the majority of the baryonic mass undergoes strong shocking so long 
as the nonlinear mass scale exceeds the Jeans mass. At infinite resolution in the Mj = 
case, it is possible that all of the baryonic material undergoes shocking. As anticipated 
from previous investigations, the highest density collapsed fraction is characteristically less 
shocked as compared with later infalling material from larger regions. The underlying cause 
for this behavior is the fact that potential wells deepen as structure grows. The highest 
density material is that which collapses earliest due to the smallest scale perturbations. 
This material falls into relatively shallow potential wells, and is only weakly shocked. 
As the structures continue to grow, progressively larger scales go nonlinear and collapse. 
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The potential wells deepen and infalling material gains more energy, resulting in stronger 
shocking and higher temperatures. 

Hydrodynamics can also play an important role in determining the distribution of the 
baryon mass, particularly in collapsed structures. In the absence of external mechanisms to 
heat the baryons (such as energy input from photoionization), during the linear phase of 
structure growth the baryons evolve as a pressureless fluid and simply follow the dominant 
dark matter. Once nonlinear collapse sets in, the baryons fall to the potential minimum, 
shock, convert their kinetic energy to thermal energy, and settle. In contrast, the dark 
matter simply passes though the potential minimum and creates a more diffuse structure 
supported by the anisotropic pressure of random velocities. This difference gives rise to 
a characteristic pattern in the baryon/dark matter ratio. Wherever the evolution is still 
linear, the baryons and dark matter simply remain at the universal mix. With the onset of 
nonlinear collapse, the baryons fall to the minimum of the potential well where they form 
a baryon enriched core, surrounded by a dark matter rich halo. We find that even in the 
absence of radiative cooling the cores of collapsed structures can become baryon enriched 
by factors of nbary/^dm ~ 2 or more, though this value is likely resolution and dimension 
dependent. If the thermal energy of the baryons is raised to the point that it rivals the 
potential energy during the collapse, the baryons will become pressure supported and stop 
collapsing at that point. In all cases we find that the dark matter is relatively unaffected by 
the baryon distribution. This is due to the fact that the dark matter dominates the mass 
density, and therefore the gravitational potential. In general it appears that under a dark 
matter dominated scenario hydrodynamics can substantially alter the characteristics of the 
baryonic material (and therefore the visible universe), such that it does not directly follow 
the true mass distribution which is dominated by the the dark matter. 

We would like to thank David Weinberg for inspiring this project, and for many useful 
discussions during its course as well. We would also like to thank the members of Ohio 
State's Astronomy Department for the use of their workstations both for performance 
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Fig. 1. — The ratio of the Jeans mass to the resolved mass (Mj/Mr) as a function of 
expansion for each of the three resolutions used in this paper (N = 64 2 , 128 2 , 256 2 ) for the 
Mj > case. The Jeans mass is calculated using the average background density of the 
universe at each expansion. 

Fig. 2. — Dark matter overdensities (pdm/pdm) for Mj = simulations. Panels arranged with 
increasing resolution along rows (N = 64 2 , 128 2 ,256 2 ), and increasing cutoff in initial input 
perturbation spectrum down columns (k c = 32, 64, 128). Part a) shows results using the full 
resolution of each simulation, while b) is calculated after resampling the simulations down to 
equivalent N = 64 2 resolution. All simulations are shown at the final time slice (expansion 
factor a/cii = 60), with grey scale intensity scaled logarithmically with dark matter density. 

Fig. 3. — Normalized dark matter overdensity distribution functions f(pdm/pdm) for Mj = 
(solid lines) and Mj > (dotted lines) simulations. Part a) shows results for full simulations, 
while part b) shows all simulations degraded to equivalent N = 64 2 resolution. 

Fig. 4. — Baryon overdensities pbary/pbary for a) Mj = at a/ai = 30, b) Mj = at 
a/cii = 60, c) Mj > at a/aj = 30, and d) Mj > at a/flj = 60. Panels arranged as in 
Figure @. 

Fig. 5. — Baryon overdensities for k c = 32 simulations, resampled to N = 64 2 resolution 
as in Figure |2|b. Note these panels represent the same simulations as the top rows of the 
previous figures. Shown are expansion factors a) a/ai = 30 and b) a/ cii = 60. Panels 
are arranged with increasing simulation resolution (N = 64 2 , 128 2 , 256 2 ) along rows, and 
increasing baryon Jeans mass (Mj = 0, Mj > 0) down columns. 

Fig. 6. — Normalized baryon overdensity distribution functions f(pbar y /Pbar y ) for Mj = 
(solid lines) and Mj > (dotted lines). We show the full resolution results for a) a/ai = 30 
and b) a/ai = 60, as well as results when each simulation is degraded to iV = 64 2 resolution 
at c) a/ai — 30 and d) a/a,, = 60. Panels arranged as in Figure [3| 

Fig. 7. — Kolmogorov-Smirnov statistic D comparing /(pbary) between simulations. Each 
line type corresponds to one simulation which is compared to each simulation listed on the 
ordinate axis, where the simulations are denoted as iV : k c . Note that the K-S statistic for 
comparing an individual simulation to itself is formally D = 0, but for the sake of clarity we 
have interpolated over these points in this plot. The panels are arranged with Jeans mass 
Mj increasing along rows, and expansion a/ai increasing down columns. 

Fig. 8. — Baryon mass distribution for the Mj = simulations as a function of overdensity 
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and temperature /(pbary/Pbary, T) (upper row) and /(pbary/Pbary, ^/T ad ) (lower row). T ad = 
T (p/po) 7_1 is defined as the temperature the gas would have due solely to adiabatic 
processes. Panels arranged as in Figure ^|. 

Fig. 9. — Average baryon to dark matter mixture as a function of baryon overdensity at 
a/di = 60. The baryon to dark matter mixture is defined as nbary/^dm = ^dmPbary/^baryPdm, 
so that nbary/^dm > 1 represents baryon enriched material, while ribary/^dm < 1 is baryon 
depleted. In each panel the solid line shows the measured average baryon to dark matter 
mixture, while the dashed lines represent the mixtures such that 10% of the mass at each 
overdensity is above and below the enclosed region. The dotted line shows the universal 
average ribary/wdm = 1- The top and bottom rows represent the Mj = and Mj > 
simulations, respectively. 
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